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Foreword 

This report summarizes the research accomplishments performed under the NASA 
langley Research Center Grant No. NAG 1-1749, entitled: “Application of the Spec- 
tral Element Method to Interior Noise Problems,” for the period August 1, 1995 to 
July 31, 1998. The primary effort of this research project was focused the develop- 
ment of analytical methods for the accurate prediction of structural acoustic noise 
and response. Of particular interest was the development of curved frame and shell 
spectral elements for the efficient computational of structural response and of schemes 
to match this to the surrounding fluid. 
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Introduction 

This report summarizes research to develop a capability for analysis of interior 
noise in enclosed structures when acoustically excited by an external random source. 
Of particular interest was the application to the study of noise and vibration trans- 
mission in thin-walled structures as typified by aircraft fuselages. 

The basic idea of the research is to reformulate the structure-fluid interaction 
problem using a matrix methodology based on the spectral element method. In 
this way, a wave analysis of the problem is retained, yet complex structures can be 
handled in a convenient manner. The analysis, which is formulated in the frequency 
domain, is capable of providing detailed information on the response (in either the 
time or frequency domains) to broad band excitation in any frequency range. The 
two significant features of the problem, namely; that the loadings on the structure 
are distributed and that the structures themselves form enclosed cavities (or cabins), 
are handled well by this formulation. 

The spectral element method is a powerful tool for wave propagation problems. It 
is a matrix method based on wave solutions that exactly satisfy the governing equa- 
tions and the boundary conditions. That is, the method exactly models the inertia 
properties and therefore the elements can be large, in fact spanning the region between 
discontinuities. Consequently, the system size is much smaller compared to a conven- 
tional element formulation. Furthermore, the effects of damping, viscoelasticity, and 
higher order structural models can be easily incorporated in the formulation. 

This report focuses on three related topics. The first concerns the development of a 
curved frame spectral element, the second shows how the spectral element method for 
wave propagation in folded plate structures is extended to problems involving curved 
segmented plates. These are of significance because by combining these curved spec- 
tral elements with previously presented flat spectral elements, the dynamic response 
of geometrically complex structures can be determined. The third topic shows how 
spectral elements, which incorporate the effect of fluid loading on the structure, are 
developed for analyzing acoustic radiation from dynamically loaded extended plates. 




Chapter 1 

Deep Curved Beams and Rings 


Using the curved beam equivalent of Timoshenko beam theory, the spectral element 
method is extended to problems involving curved members. By combining these 
curved beam spectral elements with previously presented straight spectral elements, 
the dynamic response of geometrically complex structures can be determined. Of 
particular interest in this paper is the coupling that naturally occurs between the 
axial and transverse degrees of freedom and how it affects the element formulation. 
As an example of the utility of this element, the point excitation of an infinite curved 
beam and a closed ring is demonstrated. 
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1.1 Introduction 

There is considerable intrinsic interest in waves in curved beams because of their use 
as arches, helical springs and rings, in such structures as aircraft fuselages and ship 
hulls. Some idea of the range of applications can be found in References [1, 2, 3, 4]. 
There is also the special case of negligible bending stiffness which corresponds to waves 
in cables and power lines; an interesting analysis of this is given in Reference [5] . This 
paper is a continuation of References [6, 7, 8] which developed a matrix methodology 
for analyzing wave propagation in complex frame structures. Specifically, we extend 
the spectral element method to include deep curved beam elements. 

Two elements are derived: a semi-infinite element, termed a throw-off element, 
and a finite length element, termed a two-noded element. The throw-off element is 
important in wave propagation problems since it is used to model remote boundaries 
which do not reflect waves. Both of these elements exactly model the distribution of 
mass and rotational inertia and thus can be of any length. While it is possible to 
model a curved beam as a collection of straight or curved segments, as in conventional 
element formulations [9, 10], the fact that spectral elements can be very long dictates 
that we use only one element between any two joints or points of discontinuity. 

An interesting aspect of curved beams is the coupling that occurs between the 
longitudinal and flexural degrees of freedom. The coupling is interesting in the fact 
that purely axial or transverse excitations will cause both longitudinal and flexural 
responses in the curved beam. Unlike straight beam modeling where the coupling 
between the degrees of freedom occurs only at attachment nodes, the curved beam 
possesses coupling at the differential level. That ij., the longitudinal and flexural 
motion of the curved beam are coupled through the equations of motion, and results 
in a spectrum relation that is relatively complicated. Therefore a portion of this 
paper is devoted to discussing the spectrum relation in some detail. 

While curved elements can be combined with straight elements to form geometri- 
cally complex structures, we will not emphasize that aspect of their use. Rather, we 
wish to focus specifically on some of the wave propagation aspects. We look at two 
problems: an infinite curved beam and a closed ring. The infinite beam is used to 
demonstrate the coupling between the longitudinal and flexural degrees of freedom. 
The ring illustrates how the point excitation of a simple structure can be viewed as 
either a wave propagation problem, a vibrations problem, or a rigid body motion 
problem. 
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1.2 Spectral Analysis of a Deep Curved Beam 


Consider the curved beam segment shown in Figure 1.1. Following Reference [10], 
the 2-D deformation of the beam can be approximated as 

u(s, y, t ) « it (a, t ) - y<}>(s, t ) , v(s, y, t) « u(s, t ) 


where u(s, t ) is the mean mid-plane circumferential displacement, i/(s, t) is the radial 
displacement, and <p(s , t) is a rotation about the mid-plane. This deformation leads 
to the non-zero strains 


du v d<f> _ udv 
ess= ds~R~ V !te ’ %y = R + ds ~ ^ 

Other curved beam theories have slightly different expressions for these strains; the 
present theory is closest to that of the Timoshenko straight beam [11]. The most 
significant aspect of this strain-displacement relation is the non-zero centroidal strain 
(at y = 0) even if v is the only deflection. This will give rise to the coupling of the 
two displacements. 

The strain energy for the small segment of curved beam in plane stress is 


U = i JjEil + GtIJ di = i I l [EAC 


,du v , 2 _ . r , ,ii dv ,, 2 . , 

y s ~R ] +EH Ts ) +GAK ' { R + Ts~ ,t ’ ) ]is 
where E is the Young’s modulus, G is the shear modulus; EA, GAK\ and El are the 
extensional, shear and bending stiffnesses, respectively; V is the volume and L the 
segment length. Note that, as is commonly done, we have associated an adjustable K\ 
with the shear stiffness. Typically, this it can be taken close to unity; Reference [11] 
gives a discussion of the choice of K\ for the flat Timoshenko beam. The total kinetic 
energy is 


T = \ j^p[u(s, y,t ) -I- v 2 (s , y,t)] dV = | Jjp/ In 2 + pIK 2 4> 2 + pAv 2 ] ds 

where p is the material density. Here too we introduce an adjustable parameter K 2 
to be associated with the rotary inertia. An application of Hamilton’s principle [12] 
using the variations with respect to Su, 5v, and S(j>, leads to the three governing 
equations of motion (for R=constant) 


EA 


EA 


d 2 u 
d? 
1 du 
R ds 


1 dv 
Rds 

v 

R? 

El 


- GAKi 


u 1 dv 4> 
~R? + Rds ~ R 


+ GAKi 

d 2 (f) 


1 du d 2 v d<f> 

R ds + ds 2 ds 


ds 2 


+ GAKi 


u dv 
— -j- — — d) 
R ds * 


d 2 u du 

pA W + vA m 

d 2 v dv 

pA w + ” A &t 

TJ . d 2 (j) rr^94> 

pII<2 ~di? +ri <2 ~dt ^ 
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where we have added some viscous damping rj. The associated natural boundary 
conditions are given in terms of the resultant forces 


F = 


= L a - 


4 A = EA 


du 

ds 


V = 


= f o sy dA = GAKi 
J A 


udv 
R + Ts~' t ‘ 


( 1 . 2 ) 


and moment 


M = - o ss y dA = El 

Ja 

acting on the cross-section, where the integration is over the cross-sectional area A. 
When R becomes very large, the straight deep beam and elementary rod theories are 
recovered. 

Spectral analysis assumes solutions of the form 


d<j) 

ds 


(1.3) 


«(«, *) = 53 “( s ’ u ) elWt ’ v i s 4) = , (f)(s, t) = ^2 4>(s, w)e iu ' t 

(1.4) 

where the summations are over frequency. When these are substituted into the govern- 
ing differential equations, we get a set of ordinary differential equations with constant 
coefficients. These have solutions of the form 

u(s, l j) = u 0 e ~ lks , v(s, u) = v Q e ~ lks , <^(s, u) = <f) 0 e~ iks 


where k = k{uj) is the wavenumber. In this representation, the amplitudes u Q , v 0 , 
(j> 0 , and the wavenumber are as yet undetermined. On substitution these lead to the 
homogeneous system of equations 


' Qq - GAKi/R 2 

{EA + GAKi)ik/R 

GAKi/R ' 


u 0 ' 

-{EA + GAK\)ik/R 

a 2 - EA/R 2 

ikGAKi 

i 

V 0 > 

GAK4R 

-ikGAIU 

a 3 - GAKi 


{ <t> o J 


with 

= -EAk 2 -I- pAu 2 , a 2 = -GAKik 2 + pAQ - 1 , a 3 = - Elk 2 + pIK 2 u 2 


and O 2 =uj 2 - iup/p. For a non-trivial solution, the determinant must be zero and 
this allows us to determine k. This has six solutions in all, but since only k 2 terms 
appear, there are three basic modes appearing as ±k m pairs. 

For each wavenumber k m , Equation (2.7) gives u:; the relation among the ampli- 
tudes. This is homogeneous and therefore, at best, W 2 can only get amplitude ratios. 
For example, we can solve for the remaining two terms as a function of u Q . Let us 
write the solutions at a particular wavenumber k m as 


r u 0 ' 


1 1 

Vo 

> = i 


<Po J 

m 

l ** J 
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where the symbol $ indicates an amplitude ratio. Although the vector {$} shown is 
normalized with respect to u„, it is possible for other modal vectors to be normalized 
differently. This must be done for each mode k m and hence there are six vectors. We 
choose to represent these as 


[*] = [[*a],Ms 



[ 1 1 


f 1 

f *u } 

r 1 1 


r *« 1 


r 1 


< 

1 $ 1 
1 v 1 

► i 

1 

i 1 

1 1 
1 v i 

► < 

1 

► < 

1 

> 

_ 

1 J 

1 

l J 2 

l % ) 3 

i % J 

4 

1 *+ J 

5 

l n J 

6- 


where the [ 3 x 3 ] partitions and [$#] are evaluated at +k m and —k m , respectively. 
The matrices [$4] and [<h#] are referred to as modal matrices. They are fully popu- 
lated matrices and typically are not symmetric. The normalizations are arranged so 
that if we set d> 12 = $13 = $21 = $31 = 0, the uncoupled straight beam solutions are 
recovered. 

For each mode, the corresponding amplitude u om is undetermined; to make the 
notation resemble what we have already used, we will label each of the these A, B, • • •. 
The solution for the displacements can then be expressed as 


u(s) = A$ n e- iklS + B$i 2 e“ ifc2S 

v(s) = A<S>2ie~ iklS + B<S>22e~ ik2S 

4 >(s) = A<l> 3 l e- iklS + B<S> 32 e~ ik2S 


+ • 

• • + E<h 15 e +ifc2S + F$i 6 e +ifc3S 


+ ■ 

■ ■ + E$ 25 e +ifc2S + F$ 26 e +ifc3S 


+ • 

• • + E<$> 35 e +ik2S + F$ 36 e +ifc3S 

( 1 . 6 ) 


where the terms <&ij are the amplitude ratios. The coefficients A, B, • • •, F are to be 
determined from the boundary conditions. 

It is apparent that the spectrum relation plays a central role is the solution, and 
since the characteristic equation is rather complicated, we look at its solution in 
greater detail next. 
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1.3 Discussion of the Spectrum Relations 

The characteristic equation to determine the wavenumber k is formed by setting 
the determinant of system (2.7) to zero. To simplify the expressions, first introduce 
the wavenumbers kp = pAui 2 /EA, k 2 s = pAu) 2 /GAK\, k 2 = plQ 2 / El, and k B = 
pAu 2 /EI. On expansion, we find the characteristic equation can be rearranged as 

k ® + 02 ^ T Oi k 2 + do = 0 (1.7) 


where 

0,2 = — {k"p + A ; 2 + k 2 + 2 /i? 2 ) 

a! = kpkg + k 2 P k) + kgk 2 - k B - (k 2 P + k^ - 2k]) /R 2 + l/R 4 
a 0 = (—k^ + l/R 2 )(k 2 k] — k‘ i B — k 2 /R 2 ) 

This has six solutions in all, but since only k 2 terms appear, there are three basic 
modes appearing as ±kp one associated with the longitudinal behavior and two 
associated with the flexural behavior. This can be seen by noting that for very large 
R the characteristic equation can be factored into 

(. k 2 — k 2 P )[k 4 — (kg + k 2 )k 2 — (k B - k^k 2 )] = 0 

where the term in parenthesis is the characteristic equation for the longitudinal motion 
in a rod [11] and the term in the square brackets is the characteristic equation for 
the uncoupled flexural response of a Timoshenko beam [6]. In general, of course, the 
modes are coupled and it is not proper to speak of a longitudinal mode or a flexural 
mode. 

Before we solve for the spectrum relations, it is beneficial to check certain features. 
First, we see if there is a cut-off frequency; set k = O n the characteristic equation to 
get 

a 0 = (-kp + 1 /R 2 )(k 2 s k 2 - k% - kj/R 2 ) — 0 

After setting the damping to zero, this yields the two cut-off frequencies, 

_ 1 [EA _c 0 _ [GAR] GAK X c 0 

Wcl R\j pA R ’ Uc2 V Pi + pAR 2 ~ 2 h 

where h is the depth of the beam. The presence of a cut-off frequency is typical of 
elastically coupled systems. It is interesting to note that u cl depends on the radius of 
curvature, while u c2 is dominated by the beam depth. It is clear from the expression 
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for uj c2 that for slender beams, this cut-off frequency is very large; this cut-off is 
associated with the Timoshenko second mode. 

Now look at when the frequency is zero, the characteristic equation can be factored 
as 

W - ±)(k 2 - jg) - 0 

Only one root goes through zero, the others are a double root on the real axis. 

We must solve a cubic equation in order to get the full behavior of the spectrum 
relations. The formulas for doing this are more complicated than for the quadratic 
equation and can be found in Appendix A. Figure 2.2 shows the first three spectrum 
relations, these correspond to propagating waves, and are characterized by a negative- 
only imaginary component. To emphasize the coupling, the plot is for an aluminum 
beam that is 100 mm (4.0 in) deep with a radius of curvature, R = 100mm (4.0 in). 
We clearly see the cut-off frequencies in the first and third modes. Note that k 3 has 
a negative real component at low frequencies and it might therefore be thought that 
this violates the radiation condition for waves propagation in the positive direction. 
In our approach, the wavenumbers always have an imaginary component — even 
predominantly real-only modes such as k 2 have an imaginary component arising from 
the damping. Thus the criterion is based on dissipation of energy in the positive 
direction. A negative real component is expected to lead to a standing wave. 

For the later examples, we will use an aluminum beam that is 25.4 mm (1.0 in) 
deep with a radius of curvature, R = 254 mm (10.0 in). The spectrum relations 
for this beam are plotted in Figure 1.3. Only the cut-off frequency in the first mode 
appears in the frequency range of interest. In the low frequency region (0-500 Hz), we 
see some coupling effects due to the thickness of the beam. This coupling diminishes 
as the radius of curvature increases. 

In order to better understand what the spectrum relations are doing, it is worth 
while to consider a set of approximate spectrum relations. The three roots of the 
characteristic equation with relatively large R/h ratio can be approximated nicely 
from the straight beam case as 
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where ki is the longitudinal dominated mode and k 2 and A; 3 are the flexural dominated 
modes. These approximations, also plotted in Figure 1.3, allow us to make a few 
statements about the behavior of the coupling. In comparison to the uncoupled 
modes, the major effect is in the longitudinally dominated mode. It can be seen that 
the behavior is similar to a rod with elastic constraint [11], that is, it is imaginary- 
only up to the cut-off frequency and then becomes real-only. Thus the low frequency 
components evanesce indicating a transfer of energy to the other modes. The 1/2 ft 2 
terms in k 2 and k 3 are acting like a compressive pre-stress on the beam [11], This 
has the effect that for the propagating flexural mode the group speed is increased 
indicating an effective increase in stiffness. However, the coupling generally will not 
cause drastic changes in the spectrums of the flexural dominated modes. Indeed, if 
h « .ft/100, it is clear from Figure 1.3 that the only difference is in the first mode. 
Finally, notice how the spectra are almost identical to the uncoupled spectra at the 
higher frequencies. 
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1.4 Point Excitation of a Curved Beam 

As a prelude to considering curved beams of finite length, we begin by looking at 
the point excitation of an infinite curved beam; physically this would mean that the 
beam is in the form of a helix. We will use two types of force histories, one that is 
relatively broad-banded in frequency and two that are relatively narrow-banded, in 
order to demonstrate how the coupling of the modes changes with frequency. These 
are shown, along with their normalized amplitude spectrums, in Figure 1.4. 

Our approach to the solution parallels the problems presented in References [6, 8, 
11]; what makes this case interesting is that we now have three coupled modes. The 
solution for the forward propagating terms is written as 

u(s) = A4>ne- ifclS + B4> 12 e- ifc2S + C$ 13 e^ 3S 
v{s) = A* 21 e" ifclS + B* 22 e~ ifc2S + 

j>(s) = A$3ie" lfclS + B$32e- ifc2S + C$33e^ (1.9) 

From the free body diagram of the excitation region of the infinite beam it can be 
shown that at the loading site s = 0, 

it(0, t) = 0, 0(0, t) = 0, V (0, t) = — \P{t) 

where we are considering a transverse impact. The first two of these allow the solution 
to be written as 


u(s) = A[$ne _ifclS + ik2S + /3*i3e _ifc3JS ] 

v{s ) = A[$ 21 e'^ lS + a* 22 e“ ifc2S + p<f> 23 e~ ik3S ] 
0(s) = A[<h 31 e- lfclS + a $32e“ ife2S +/3$33e- ifc3S ] 


where the coefficients a and j3 are given by 

*11*33 — *31*13 


a = — 


*12*33 — *32*13 
We determine A from the shear relation 

_ip = V = GAKi 

After differentiation, this leads to 

-P 


/? = 


* 11*32 — * 3 . *12 
* 13*32 - * 33 * 


33^12 


u dv 

n + lTs~ 


A = 


2GAK\ q\ + aq2 + f3q^\ 


Qm — ^1^2m ^3m 
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The solution is arranged so that if we set R = oo then $12 = $13 = <f> 21 = $31 = 0, 
and the uncoupled straight beam solution is recovered. Care must be taken, however, 
in order to approach the proper limit since both a and beta approach infinity. It turn 
out that a/(3 = 1 in the limit and we get, for example, 


<i>(s) = 


-p 


2 Em - ki) 


e 


~ik2S 




Note that the spectrum relations also change. 

The velocity reconstructions for the broad band input are shown in Figure 3.2. 
There are two points of interest. First, note how the initial zero axial velocity even- 
tually becomes significant. Second, note the oscillatory behavior of the transverse 
velocity. Although the force excitation lasts only about 200 fxs, the beam near s = 0 
continues to oscillate in an almost resonant like fashion. Actually, a standing wave 
has been established. The figure also shows the separate contributions from each 
mode for the response v{s = 3 R, t). It is clear that the first mode is contributing the 
ringing behavior. 

This is more evident when we look at these velocity responses in the frequency 
domain. It is clear in Figure 2.8 that there is a peak in the u-velocity response which 
corresponds with the cut-off frequency. Furthermore, this peak is coming entirely 
from the first mode. The broadband excitation has identified the cut-off frequency in 
the first mode. 

A final point of interest for this example looks further at the effect of the cut-off 
frequency in the longitudinal mode. Below this frequency there is only one propa- 
gating mode in the beam, while above the frequency two propagating modes exist. 
To illustrate this point, we excite the curved beam with the narrow banded force 
histories shown in Figure 1.4. The pulses are chosen so that they just bracket the 
cut-off frequency. Figure 1.7 clearly shows the presence of the second propagating 
mode above the cut-off frequency. Note that the two propagating modes are present 
in both the transverse and longitudinal responses b it the amplitude ratios are dif- 
ferent. This plot is emphasizing the nature of the solution of Equation (2.8) as a 
collection of mode responses. 
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1.5 Spectral Element Formulation 

As seen from the point excitation example of the last section, the formulation, al- 
though relatively straight forward, requires a significant amount of manipulation. A 
similar approach for connected curved beams or even for a finite curved beam would 
be very cumbersome, hence we now develop a matrix formulation to facilitate these 
manipulations. 

Consider a segment of curved beam of length L. We begin by expressing the 
displacements as 

u{s) = A<S> n e- iklS + B<f> l2 e~ ik2S + --- + E$ l5 e- ik2 ( L ~ s) + F<Z> l6 e- ik3iL ~ s) 

v(s) = A<f> 2l e- ikxs + B<f> 22 e~ ik2S + --- + E$2se- ik2{L - s) + F<f>26e~ ik3{L ~ s) 

i(s) = A$3ie“ ifclS + B$ 32 e“ ifc2S + • • • + E$ 35 e" ifcl(L- ' ) + F^ee^^l-lO) 

The length is introduced to include reflections coming from a boundary located at 

s = L. This displacement solution can now be re-written as 

{«(*), »(«), M) t = {w}W 

= {4>}iAe-'*‘ s + ■ • ■ + {4} 6 Fe- itl(i - s) 

We will re-write this in an even more compact matrix form; so as to make the matrix 
notation a little more accessible, we will take the developments of the rod as the 
archetype and use its notation (except changed to matrices). The 1-D solution for a 
rod [11] is represented as 

u(x) = [e-'* lI ]A + [ e - ifcl(L - x) ]B 

where A and B are associated with the forward moving and backward moving waves, 
respectively. The displacement for the curved beam is written as 

{«}« = l^ir«(*)J{A} + Nfe(i - s)J{B} 

were [$4] or [<&b] are the [3 x 3] partitions of [ $ ] and 



r e- ikis 

0 

0 ■ 




r d i 

\e(s)\ = 

1 

0 0 

g-ifc 2 s 

0 

0 

e ~ik3S 

, {A} = < 

Is!' 

{B} = < 

E 

l F J 


It is the presence of the amplitude ratios that is the most significant difference. 

We wish to replace the vectors {A} and {B} in terms of the nodal displacements 
at s = 0 and s = L. That is, we introduce 

u(0) = tii , t>(0) = i>i , 0(0) = 0i ; 


u (L) = u 2 , v(L) = v 2 , 0(L) = 0 2 
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We write this in matrix notation as 


{«!, v u 4 >x) t = {u}i = {w}(s = 0) 
{u 2 , v 2 , ^ 2 } T = W 2 = W( s = L) 


= M[e(0)J{A}-[^ s ]re(L)j{B} 
= M[e(L)J{A}-Mre(0)j{B} 


Let us write all six equations as 



Solving for the coefficients gives 


/ {A} 1 _ r n — 11 ( {“}i 1 _ r r 1 f { u }i 1 _ [[G 11 ] [G 42 ]] f {^}i \ 

\{B}/" iy J l{u} 2 /- lGJ l{u}2/ _ [[G , 21 ] [G 22 }\ l {«} 2 J 


where each partition of [ G ] is of size [3x3]. We are now in a position to write the 
displacements in terms of the shape functions. They are 


(W}(s) = b(s)]i{“}i + (1.11) 

where the [3 x 3] matrix of shape functions are defined as 

b( s )]i = Mr^)J[Gn] + Mre(T-s)j[G 21 ] 

b( s )]2 = [^. 4 ] fe(s)J [G 12 ] + fe(Z/ — s)J [G^] (1.12) 


There are a total of 3x3x2 = 18 shape functions in all. While not obvious from the 
above, it turns out that, even in this general case, the collection of shape functions 
associated with the degrees of freedom at the seconl node are the mirror image of 
the first set. 

Figure 3.6 shows the g 22 (s) shape function of a 270° beam segment at a number 
of frequencies. This shape function is associated with the ?>i degree of freedom and 
thus can be plotted as a radial displacement off the original shape. It is similar to 
conventional shape functions except that it is frequency dependent, typically com- 
plex, and can represent the behavior of very large beam segments. The other shape 
functions behave in a similar manner. It is clear from this figure and from the above 
developments, that once the nodal degrees of freedom are determined that the shape 
functions can be used to compute the responses at any intermediate locations. This 
is a crucial attribute since the beam segments or elements can be very long. 

The next step in the element development is to c erive a stiffness relation for the 
beam segment. The process is simply that of expressing the resultant forces and 
moments from Equations (2.5) and (1.3) in terms of Ihe displacement solutions given 
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in Equation (2.12). We can write resultants associated with the boundary conditions 
in terms of the displacements in matrix form as 

MM = {f, v, m} t (s) = ( a |{w}(») 

where [ d ] is the matrix collection of differential operators of size [3 x 3]. After 
substituting for {£/}(s) in terms of the shape functions get 

OTM = [a]b(s)]lMl + [9][^)]2{«} 2 

= [dtfOOliMi + [d0(s)] 2 M 2 (1.13) 

Relating the member resultants at s = 0 and s — L to the nodal loads at the same 
locations leads to the stiffness relation 

/ {F)i \ = [ l-9»(0)]i [-8 j(o)] 2 1 r {tiji l 

HfW l[+8j(i)], |+8j(t)]2jU“h/ 1 ' 

or simply 

{F} = [^(w)]{u} , {F} = {Fi, Vj, Mi; F 2 , V 2 , M2 } T , {«} = {^l, Vi, <^i; U 2 , i>2, 0 2 } r 

where [ k ] is the [6 x 6] dynamic element stiffness matrix. This stiffness matrix is 
frequency dependent, complex, and symmetric. 

For illustrative purposes, Figure 2.4 shows the normalized k n , k 22 , and k 33 diago- 
nal terms for the same beam used to illustrate the shape functions. As is typical with 
spectral elements, they exhibit a very large dynamic range. The normalizations are 
with respect to the stiffnesses for straight thin beams [12]. That is, they are presented 
as 

k n /(EA/L) , k 22 /(12EI /L 3 ) , k 33 /{\2EI/L) 

Note that the k n stiffness is substantially less than the elementary rod values, and 
it is only after the cut-off frequency (« 3 kHz ) does it become greater than unity. 

This can be understood by considering the static (low-frequency) axial loading of the 
beam fixed at one end — because of the curvature, the axial load creates a moment 
and consequently, the beam has a great deal more flexibility than the corresponding 
straight case. 

The stiffness relation for the throw-off or semi-infinite element is simply 

{F}i = [-dgiOMuh = fauOKu}, , to(s)], = [^][e(s)J[G„] (1.15) 

The stiffness matrix in this case is a [3 x 3] frequency dependent, complex and sym- 
metric matrix. This can be used to recover the results presented for the impact of an 
infinite or semi-infinite curved beam 
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1.6 Point Excitation of a Closed Ring 


The advantage of the curved element is that it can be combined with other elements, 
either flat or curved to form significantly more complex structures. References [6, 13] 
discuss the programming structure required for multiply connected spectral elements. 
However, we will not emphasize this aspect of their use, rather, we will perform the 
example of the point excitation of a closed ring. The example demonstrates the basic 
method of joining elements in addition to demonstrating the utility of the spectral 
method. 

We form the closed ring by combining two curved elements of the same material 
properties each of which is a half-circle. Because of symmetry, we could model the 
ring with a single half-circle element or as a single whole-circle element. In the latter 
case half of the applied load must be placed at each node and the extra conditions of 
ui = 0, (f>i = 0; u 2 = 0, 4>2 = 0 be imposed. The two element model was chosen so as 
to verify the assemblage process. 

The elements can be joined at the common nodes by merely summing together 
the appropriate dynamic stiffness matrix components to form a global stiffness matrix 
[ K ]. We must therefore first rotate each element stiffness to this global system. The 
transformation requires the use of a simple [6 x 6] rotation matrix [ T ] that is of the 
form 


[T} = 


[*(«)] o 

0 [R{a + :>P)] 

where R(6 ) is the [3 x 3] rotation matrix [12], This ;akes into account that the ends 
of the curved element are oriented differently to each other. In the plane case, if the 
nodal coordinates are (aq, y,) and (x 2 , y 2 ), then with 2 p = yj(x 2 - aq) 2 + (y 2 - y^ 2 
and q = \J R' 2 — p 2 , the angles are given by 


a = tan’ 


( x 2 ~ x i)p — ( 2/2 - Vi )?' 


0 = tan 


i 


i x 2 - X\ )q + (V2 - 2/1 W ’ 

Since only two elements are present, we will assign the coordinate system of one to 
be the global coordinate system and rotate the other. 

The matrix [ T ] transforms the vectors of nodal displacements and nodal forces 
to the global system as follows 


{F} = [ T ]{£} , {«} = [r]{&} 


where the barred quantities represent the local coordinate system. As a consequence, 
the element stiffness matrix in global coordinates ca i be written as 


[ k } = [ T ] T [ k ][ T | 
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The global stiffness matrix is then simply formed by adding the two global element 
stiffness matrices together. 

With the elements connected, we can now determine the responses of the closed 
ring due to a point excitation. The velocity reconstructions for three points, 9 — 
0°, 90°, 180°, along the ring are shown in Figure 2.5. For comparison purposes, the 
same ring was modeled using 64 straight Timoshenko spectral elements. The two sets 
of results are indistinguishable. 

Unlike the previous example of an infinite curved beam, there are multiple re- 
flections occurring. The waves are traveling around the ring and interacting with 
each other. From Figure 2.6, we can see that many resonance frequencies are being 
established. These would be the frequencies of interest if a vibration experiment were 
being performed. It is interesting to note that the v response at 90° is missing many 
of the intermediate resonances, but these are present in the u response. Furthermore, 
we see the effect of the first mode cut-off frequency in the v responses but not in u. 
This is consistent with Figure 2.8 for the infinite beam. 

As a final illustration of the results, we integrate the velocity responses for the 
ring to determine how it is deforming over time. Figure 1.12 illustrates the exag- 
gerated deformation (multiplication by 300) of the ring due to the point excitation. 
It is interesting to note the rigid body motion that is occurring although it was not 
specifically addressed in the solution. On the time scale shown, the problem reduces 
simply to a transfer of momentum. The force transfers an impulse of about 0.31 Ns 
to the ring of mass 2.80 kg, thus the force should cause the ring to move at about 
113mm/s. From the figure, we can approximate the ring velocity as 103mm/s, the 
difference being due to the damping present in the modeling. 

The flexibility of the spectral element approach has allowed us to view the prob- 
lem of the impact of a ring as one of wave propagation, or vibration, or rigid body 
dynamics. This comes about because of the convenience in alternating between the 
time and frequency domains. 
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1.7 Discussion 

A deep curved beam element was developed that extends the variety of problems the 
spectral element method can handle. One of the challenges that arose during the 
development was the problem of determining the spectrum relations. Even though 
cubic solvers are well known, it is not always certain which branch should be chosen 
after a branch point occurs. In the quadratic case, which occurs for straight beams, 
the ambiguity is removed by adding damping. This becomes essential in the cubic 
case so that the phase of each k m {uj) can be tracked correctly. The effect of the 
damping is to separate each of the modes. With this separation, we can also reliably 
identify the associated amplitude ratios. 

The element developed possesses all of the features of the spectral element method 
which make the method desirable for solving dynamic problems. Primary among these 
are that the element can be very long and that the frequency domain formulation 
allows the system response functions to be determined automatically. This latter 
attribute, along with the fast Fourier transform, enables a duality between the time 
domain and frequency domain to be presented conveniently. Information not readily 
seen in one is often detectable in the other. The example problems demonstrate these 
features. 
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Appendix A: Computing the Spectrum Relations 

The formulas used to solve the characteristic equation are based on References [14, 
15]. By replacing k 2 with 2 , the characteristic equation can be re-written as 

2 3 -I- a? z 2 4- o,\z + a 0 = 0 

where the coefficients are insured of being complex by adding damping to the system. 
Now compute the terms 

Q — i a 2 ~ 3ai)/9 , r = (2a 3 — 9a 2 ai -I- 27a 0 )/54 

and 

si = — [r + iy/V 2 — <y 3 ] , 

The three roots are then given by 

Z\ = (si + s 2 ) - |a 2 

Z2 — — 2 ( s l + s 2) — 5^2 + * 2 ( s l — S 2 )\/3 

Z 3 = -^(s 1 + s 2 )-ja 2 -i^(s 1 -s 2 )\/3 

The six spectrum relations are given as ±\/z. 

A difficulty arises in choosing the appropriate square or cube root since we wish 

to compute single Zj at a time. The issue is that while the n th root of a complex 

number 2 is given by 

2 = a + ib = Ae 1 * , z l/n = A 

the phase = tan -1 (5/a) has an ambiguity of Nir. We remove this ambiguity by 
keeping track of the total phase. That is, starting with some value and with reference 
to the unit circle in the complex plane, as the wavenumber goes from the 4 th to the 
1 st quadrant the total phase is increased by 2ir. Conversely, if the wavenumber goes 
from the 1 st to the 4 th quadrant the total phase is decreased by 2i r. The n th root is 
then given by 

z l fn — j^ l ln e i{<t‘+NK)ln 

This scheme works quite robustly when the frequency increment is not too large. 


«2 


r + \fr 2 — <7 3 ] ^ 
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Figure 1.1: Coordinate system and resultant lozds for a deep curved beam. 
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Figure 1.2: Spectrum relations kh for a curved beam with R = 100 mm and h = 
100 mm. 
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Figure 1.3: Spectrum relations kh for a curved be<.m with R = 254 mm and h = 
25.4 mm. 
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Figure 1.10: Velocity reconstructions for the point excitation of a closed ring. 
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Figure 1.11: System response functions for the point excitation of a closed ring. 
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Figure 1.12: Deformed shapes (displacements scaled by 300) of a ring due to a point 
excitation. Time intervals are every 750 fis. 



Chapter 2 

Long Segmented Cylindrical Shells 


The spectral element method for wave propagation in folded plate structures is ex- 
tended to problems involving curved members. By combining these curved spectral 
elements with previously presented flat spectral elements, the dynamic response of 
geometrically complex structures can be determined. Of particular interest in this 
paper is the coupling that occurs naturally between the in-plane and transverse de- 
grees of freedom and how it affects the element formulation. As an example of the 
utility of this element, the point excitation of an infinite curved shell and a closed 
cylinder is demonstrated. 
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2.1 Introduction 

There is considerable intrinsic interest in waves in curved members because of their 
use in such structures as arches, containment vessels, aircraft fuselages, and ship hulls. 
Some idea of the range of applications is described by Gould (1988) and thorough 
treatments of their formulation are given by Leissa (1973) and Markus (1988). 

A frequency domain matrix methodology for analyzing wave propagation in com- 
plex folded plate structures was developed by Rizzi and Doyle (1992), Danial and 
Doyle (1995), and Danial et al (1996). This paper is a continuation of those researches 
but now extended to include curved segmented shells. At present, we consider only 
circular uniform cylinders but the segments can be of arbitrary length in the hoop 
direction. (In this way it is different from the finite strip method described by Hinton 
et al (1995).) Two spectral elements are derived: a single noded semi-infinite throw- 
off element, and a finite length two-noded element. Both of these elements exactly 
model the distribution of mass and thus can be of any length. 

A flat plate element has eight degrees of freedom (Danial et al, 1996) and the 
assembled [8 x 8] stiffness is achieved as a combination of two [4 x 4] elements plus 
a rotation matrix. But the curved element has all eight degrees of freedom coupled 
and hence we must tackle the stiffness matrix directly as an [8 x 8] system. This is 
too cumbersome to do explicitly; consequently, we lay out a computer based method 
for establishing the shape functions and subsequently, the stiffness matrix. This 
adds to the computational burden, but as compensation we get an approach that is 
conceptually simpler and helps to unify the special results established in the earlier 
references. 

An interesting aspect of shells is the coupling that occurs between the in-plane 
and flexural degrees of freedom. Unlike folded plate modeling where the coupling 
between the degrees of freedom occurs only at attachment nodes, the curved shell 
segment possesses coupling at the level of the differential equations of motion. This 
results in a spectrum relation that is relatively complicated. Therefore a portion of 
this paper is devoted to discussing the spectrum relation in some detail. 

While curved shell segments can be combined with flat elements to form geomet- 
rically complex structures, we will not emphasize that aspect of their use. Rather, 
we wish to focus specifically on some of the wave propagation aspects. We look at 
two variations of a closed cylinder problem. The first is used to verify the accuracy 
of the formulation while the second explores the nature of the wave reflections in the 
hoop direction. 
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2.2 Spectral Analysis of a Cylindrical Shells 

There are a variety of statements of the governing equations for shells; to be consistent 
with the Danial et al (1996) paper, we briefly summarize the derivation. To this end, 
we find it most expedient to first specify the deformation, obtain the strain and kinetic 
energies, and then use Hamilton’s principle to derive the equations of motion and the 
appropriate boundary conditions. 

Consider the segment of cylindrical shell shown in Figure 2.1. The shell has a 
radius R, thickness h, and is considered long in the y-direction. The 3-D deformation 
is approximated as 

u(s,y,z) « u(s,y) - z 
- / \ . . dw 

v{s,y,z) « v{s,y)-z— (2.1) 

w(s,y,z) « w(s, y) 

where u(s, y ) is the mean mid-plane circumferential displacement, v(s, y) is the length- 
wise displacement, and w(s,y ) is the radial displacement. This deformation leads to 
the non-zero strains 


w du 

-R + lte~ Z 


d 2 w 1 du 
'ds 2 + HYs 


dv du ( 
Ts + di ~ z { 2 


d 2 w 1 du 
dsdy Rdy 


Other shell theories have slightly different expressions for these strains; the present 
theory is closest to that of Reissner (1946) and Naghdi (1964). Excellent surveys 
of the different theories are given by Leissa (1973), and Markus (1988). The theory 
developed here is the shell equivalent of the classical plate theory and the Bernoulli- 
Euler beam theory. 

An application of Hamilton’s principle (Doyle, 1991) leads to the three governing 
equations of motion (for R = constant ) 


i Wi x d 2 u _ Wi , d 2 v 1 dw, 
E ^ds 2 + 5(1 ^ dy 2 5(1 + U) lhdy ~ R'ih^ 
1 f <9 2 u j d 2 u d^w d 3 w.. 

+ R? D ^ds 2+2 ^ ~ U ^d^ + R ^d^d~ S + Ts^ 
p.i,,, N d 2 u ,d 2 v , d 2 v V dw. 

£[ 5 ( ‘ + U) d^y + W 5< V) W “ R~itf 


d 2 u du 

ph w +,,h ai 

d 2 v dv 

+ ’> h m 
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1 du v dv w . 

~ E ^Rd^ + Rdy~ ^ 

^ r d A w d*w d 4 w 1 r dPu 8Pu , d 2 w dw 

+D \-^j + 2 b + ~^rz 1 + ~s D \. ~~P h ~^~^ h ~ 


(2.3) 


ds 4 ' “ ds 2 dy 2 ' dy il ' R~ v dy 2 ds ' ds 3> r " dt 2 1 dt 

where u is Poisson’s ratio, D = Eh 3 / 12(1 — zT 2 ), and E = Eh/( 1 — z/ 2 ). We have 
added some viscous damping through the y terms. 

Let the virtual work of the applied loads be 

dw 

SV = -Q u Su - Q v 6v - Q w 5w - Q^oip , ip — — 

then Hamilton’s principle also gives the associated boundary conditions on the side 
s = constant as: 


u 


v 


w 


, dw 

'*= d7 


or 


or 


or 


or 


Qu ~ E 

Qv = i(l - u)E 


du w dv 
ds ~ R + V dy 


+ l R° 


d 2 w 1 du d 2 w 1 

ds 2 R ds dy 2 


dv du 
ds dy 


Qw = i? 

QtP — D 


-D 


d 2 u . <9 2 u 

+ (1 - ") 


ds 2 


dy 2 


- D 


dPw . cPw 

+ ( 2-z/)- 


5s 3 


dsdy 2 


d 2 w d 2 w 1 

9s 2 ^ 9z/ 2 /? 5s 


(2.4) 


We can form the resultants per unit length as 
N ss — J* dz , — J" vj dz , Hfjs — 'J' (? SS Z dz , 


M, 


SJ/ 


-/ 


<7^2,2 


where the integrations are over the cross-section. After substituting for the stresses 
and strains in terms of our approximations, we get that the natural boundary condi- 
tions are equivalent to specifying 

Qu — Nss 4” ~ft^ ss 


Qv = N, 
Qw 


sy 

dM s 


dMsy y 

6 — v sz 


(2.5) 


ds ~ dy 
Qip — lUss 

The first of these resembles the resultant load expression used for curve beams (Bilodeau 
and Doyle, 1997), while the third resembles the Kirchhoff shear stress relation (Doyle, 
1997). 

Spectral analysis assumes solutions of the form 

2nn 


N O 7 r tti 

u(s, y,t)=J2J2 u(s, u? n ,£ m ) e -' Uy e' UJ ’ lt , = 

n=l m= 0 


w 


U)n — 


(2.6) 
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where W and T are the space and time windows, respectively. The use of these 
representation is documented by Danial and Doyle (1995). In this paper, we limit 
ourselves to problems that are symmetric about y = 0 and use cos (£ m y) with u and w, 
and use sin(£ m y) with v. When these are substituted into the governing differential 
equations, we get a set of ordinary differential equations with constant coefficients. 
These have solutions of the form 

u(s,u,Z) = u 0 e~ lks , v{s,u,£) = v 0 e~ iks , w(s,uj,Z) = w 0 e~ iks 

where k — k{u,Q is the wavenumber. In this representation, the amplitudes u 0 , v a , 
w 0 , and the wavenumber are as yet undetermined. On substitution, these lead to the 
homogeneous system of equations 

’ - [A ; 2 + (1 - u)e\D/R 2 - 7 [E_+ (k 2 + ?)D]ik /< Rl [V 

7_ a 2 vE£/R <u o >=0 (2.7) 

. [E + ( k 2 + e)D]ik/R -uEZ/R a 3 + E/R 2 J 

using 

a i = ~E[k 2 + |(1 - v)?] +phu 2 , j = |(1 + u)EikZ 
a 2 = + |(1 - ^)A: 2 ] + phui 2 , c *3 = D[k 2 + ^ 2 ] 2 - phw 2 

and u 2 = u 2 — iurijp. The aj, q 2 and 7 terms alone define the flat membrane 
problem, while a 3 alone defines the flat plate flexural problem. All the other terms 
are couplings. For a non-trivial solution, the determinant must be zero and this allows 
k to be determined. This has eight solutions in all, but since only k 2 terms appear, 
there are four basic modes appearing as ±kj pairs. 

For each wavenumber k jt Equation (2.7) gives the relation among the amplitudes. 
This is homogeneous and therefore, at best, we can only get amplitude ratios. For 
example, we can solve for the remaining two terms as a function of u Q . In anticipation 
of a later need, we add ip 0 = —ikw Q and write the solutions at a particular wavenumber 
kj as 



where the symbol $ indicates an amplitude ratio. Although the vector {<£} shown is 
normalized with respect to u 0 , it is possible for other modal vectors to be normalized 
differently. This must be done for each mode kj an i hence there are eight vectors. 
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We choose to represent these as [ <J> ] = [[$4], [$#]] where 


[*a] = 



r 

1 
$ 

) 




and [4>fl] is the same but evaluated at the wavenumbers — kj. That is, the [4 x 4] 
partitions [4>^] and [4>e] are evaluated at +kj and —kj, respectively; they are fully 
populated and typically are not symmetric. The normalizations are arranged so that 
the uncoupled flat plate solutions are easily recovered. 

For each mode, the corresponding amplitude u Q j is undetermined; to make the 
notation resemble what we have already used, we will label each of the these A, B, • • •. 
The solution for the displacements can then be expressed as 


u(s) = A<t> u e~ iklS + B4> 12 e- ifc2S + • • • + G$ 17 e +i * 3S + H$ 18 e +ifc4S (2.8) 

with similar expressions for v, w, and ip, but involving the amplitude ratios d> 2 j, $3 j, 
and <h 3 j, respectively. The coefficients A,B,---,H are to be determined from the 
boundary conditions. 

It is apparent that the spectrum relation k{u, £) plays a central role in the solution, 
and since the characteristic equation is rather complicated, we look at its solution in 
greater detail next. 
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2.3 Discussion of the Spectrum Relations 

The characteristic equation to determine the wavenumber is formed by setting 

the determinant of system (2.7) to zero. On expansion, we find the characteristic 
equation can be rearranged as 

k + A 3 k 6 + A 2 k^ + Aik 2 + Aq — 0 (2-9) 

The expressions for the coefficients A n are too complicated to be listed here. This 
has eight solutions in all, but since only k 2 terms appear, there are four basic modes 
appearing as ±kf two associated with the in- plane behavior and two associated with 
the flexural behavior. In general, of course, the modes are coupled and it is not proper 
to speak of a membrane mode or a flexural mode. 

We must solve a quartic equation in order to get the full behavior of the spectrum 
relations. While the formulas for doing this are relatively straight-forward, they cause 
some difficulties which are worth discussing. 

Following Abramowitz and Stegan (1965), we first write the characteristic equation 
as the conjugate factorization for z = k 2 

[z 2 + (o + 6)z + (c + d)j ^z 2 + (a — b)z + (c — d)j = 0 

This allows z to be determined from a sequence of quadratic equations. The coeffi- 
cients appearing in the above are 

a = \A?> , b — |^/ 8 c + A 3 — 4A 2 , d = \Jc? — Ao 

and c is chosen as a solution of the cubic equation 

c 3 — 2 [•^ 2 ]^ T |[A!A 3 — 4A 0 ]c T |[AqA 2 + A 2 — 4AoA 2 ] — 0 

The formulas to solve the cubic equation are basel on those of Abramowitz and 
Stegan (1965) and Press et al (1992). We first write the equation as 

c 3 -f a 2 c 2 + a x c -f Oq = 0 

and then compute the terms 

q = (al - 3ai)/9 , r = (2a 3 - 9a 2ai + 27a 0 )/54 , s 1>2 = -[r ± yV 2 - r/ 3 ] ^ 

The three roots for c are then given by 

Ci = («i + s 2 ) — |a 2 

C 2 = -|(si +S 2 ) - |a 2 + z^(,q -s 2 )\/3 

c 3 = -J(ai +«2)-Jaa -*{(■’! -*2)V3 
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The eight spectrum relations are given as ±\fz. 

These equations can be easily programmed. The difficulty that arises is in choosing 
the appropriate square or cube root in both the cubic and quadratic solvers, as well 
as the appropriate root from the cubic solver. It is seen that while there are only four 
z roots, the formulas give a multiplicity of 144 roots. This is especially acute since 
we wish to compute a single Zj at a time. The issue is that in computing the n th root 
of a complex number z = a + ib — Ae*^, the phase (j) = tan _1 (5/a) has an ambiguity 
of Ntt. 

We remove this ambiguity by keeping track of the total phase. That is, starting 
with some value and with reference to the unit circle in the complex plane, as the 
wavenumber goes from the 4 th to the 1 st quadrant the total phase is increased by 2ir. 
Algorithmically, we compute complete modes separately at each m. Starting at a large 
value of frequency, we work toward the lower frequencies, keeping track of the phases. 
The root c is chosen as the one that had the largest real value initially. Approximate 
spectrum relations are useful here in identifying the appropriate modes. This scheme 
works quite robustly when the frequency decrement is not too large. Periodic checks 
with an iterative root finder helps confirm the correctness of the roots. It must be 
said that determining the spectrum relations is now a significant operation in itself, 
and unlike the previous reported cases, they are no-longer computed on-the-fly as 
part of the structural analysis. 

An idea of the variety of behaviors is shown in Figure 2.2 for a value of £ = 2irm/W 
with m = 40. In this and the remaining examples, we consider an aluminum shell 
that is 25.4mm (1.0m) thick with a radius of curvature, R = 254mm (10.0m). Also, 
we take W = 20 m (800 m). Using these values, the curvature has a significant effect 
on the spectrum relations. The figure shows for the first four spectrum relations; 
these correspond to propagating waves, and are characterized by a negative imaginary 
component. 

Not surprising, there are many branch points in the spectrum relations. We clearly 
see three cut-off frequencies but what is interesting is that two of them are associated 
with k 2 , the in-plane shear dominated mode. Note that both k 2 and k 4 have a 
negative real component at low frequencies and it might therefore be thought that 
this violates the radiation condition for waves propagation in the positive direction. 
In our approach, the wavenumbers always have an imaginary component — even 
predominantly real-only modes such as k 3 have an imaginary component arising from 
the damping. Thus the identification criterion is based on dissipation of energy in 
the positive direction. A negative real component is expected to lead to a standing 
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wave. 

It is difficult to get a clear idea of the complete (i.e., as m varies) spectrum relations 
from Figure 2.2, it is useful, therefore, to introduce an approximation that will help 
to delineate the separate contributions. An approximation that is reasonable for thin 
shells of large radius is given by 


i.2 _ U 2 c2 

K l — K P ~ R2 S 

L2 _ 1.2 _ c2 

™2 "S S 


kl 




2R? 

3 

2R 2 


+ \j k % + AR4 R2 


'ki + 


4i? 4 


il 

R 2 


( 2 . 10 ) 


where kj, = u 2 p(l - v 2 )/E, k 2 s = u 2 p2(l + v) / E , and k 4 B = u 2 ph/D. The first 
two are the membrane dominated modes and recover the plots shown in Rizzi and 
Doyle (1992) when R is very large. The third and fourth modes are the flexural 
dominated modes and for large R they recover the flat plate spectrum relations given 
by Doyle (1997). In comparison to the uncoupled modes, the major effect of the 
curvature is in the longitudinally dominated modes because the cut-off frequencies 
are affected and hence delay the formation of real-only wavenumbers. 

Since the spectrum relations are relatively insensitive to R at large frequency, 
Then these approximations are quite useful in identifying the modes as determined 
by the root solver. 
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2.4 Spectral Element Formulation 

When the number of degrees of freedom is large, we need an organized way to handle 
establishing the shape functions and the element stiffness. The cylindrically curved 
shell segment is such a case and we take this opportunity to develop an appropriate 
matrix scheme. 

Consider a segment of shell of length L in the hoop direction. We begin by 
expressing the displacements as 

«(s) = A<f> u e~ ifclS + B<f> 12 e-* 2S + • • • + GS^e"* 3 ^" 4 * + 

with similar expressions for v , w, and ?/>. The length is introduced to include reflections 
coming from a boundary located at s = L. This displacement solution can be re- 
written as 

{u(s), u(s), w(s), i>{s)} T = (W}(s) 

= {$>! Ae~ iklS + • • • + {$} s He~ iki{L ~ s) 

We will re-write this in an even more compact matrix form; so as to make the matrix 
notation more accessible, we will take the developments of the rod as the archetype 
and use its notation (except changed to matrices). The 1-D solution for a rod is 
represented by Doyle (1997) as 

u(x) = [e~ ifclI ]A+ [ e - ifcl(L “ x) ]B 

where A and B are associated with the forward moving and backward moving waves, 
respectively. The displacement for the shell segment is written as 

{«}(») = Mr«(*)J{A} + MML - s)J{B} 


were [$^4] or [4 >b] are the [4 x 4] partitions of [ 4> ] and 
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e ~ik 4 s 


. D . 


. H , 


We wish to replace the vectors {A} and {B} in terms of the nodal displacements at 
s = 0 and s = L. That is, we introduce 


fi(0) = u x , v(0) = t>i , u)(0) = w x , V>( 0) = V’l 
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with similar, but subscripted ‘2’, terms at s = L. We write this in matrix notation as 


{«ii fiii fib, <M T = {£}i = {g}($ = o) 
{u 2 , v 2 , w u 0 2 } T = {«}2 = {U}{s = L) 


= [*Ai re(0)J {A} - [* B ] \e(L ) J {B} 
= Nre(i)J{A}-NFe(0)J{B} 


Solving for the coefficients gives 



f[Gn] 

.[G 2i ] 


[^12]] f {“}i \ 

[G 22 ] . 1 (fi}2 / 


where each partition of [ G ] is of size [4x4]. We are now in a position to write the 
displacements in terms of the shape functions. They are 


{W}(s) = [ff(«)]i{fi}i + b(s)Mfi}2 (2.12) 

where the [4 x 4] matrix of shape functions are defined as 

b(s)]i = [^A]r e ( s )J[G'n] + M[e(L-s)J[G 21 ] 

b( s )l 2 = [ < J > A]fe(s)J[G’i 2 ] + [$a][e(L - (2.13) 

There are a total of 4 x 4 x 2 = 32 shape functions. While not obvious from the 
above, it turns out that, even in this general case, the collection of shape functions 
associated with the degrees of freedom at the second node are the mirror image of 
those associated with the degrees of freedom at the first node. The results for the 
throw-off element are simply are simply those associated with [Gn]. These formulas 
can be used to recover all the shape functions already derived in the cited references. 

By way of example, Figure 3.6 shows the g 33 (s) cos(£ m y) shape function of a 270° 
shell segment; this shape function is associated with the w\ degree of freedom and 
thus can be plotted as a radial displacement off the original shape. These shape 
functions are similar in concept to conventional finite element shape functions except 
that they are frequency and wavenumber dependent, are typically complex, and can 
represent the behavior of very large segments. The other shape functions behave 
in a similar manner. It is clear from this figure anc from the above developments, 
that once the nodal degrees of freedom are determined then the shape functions can 
be used to compute the responses at any intermediate locations. This is a crucial 
attribute since the segments can be very large. 

The next step in the element development is to derive a stiffness relation for 
the shell segment. The process is simply that of expressing the resultant forces and 
moments from Equation (2.5) in terms of the displacement solutions given in Equa- 
tion (2.12). We write these resultants in matrix form as 

(^}(s) = {N ss + |m ss , N sy , V S2 , M ss } r (s) = [ d }{U}(s) 
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where [ d ] is the matrix collection of differential operators of size [4 x 4], After 
substituting for {W}(s) in terms of the shape functions get 

WOO = [ d Ms)]iMi + [ d HsOOMuh 

= [ds(s)]i{^}i + [ac?(s)] 2 {u} 2 (2.14) 


Relating the member resultants at s = 0 and s — L to the nodal loads at the same 
locations leads to the stiffness relation 

{F},1 = [[-35(0)1, [-95(0)k 

{F)J L[+9s(F)], [+9s(Ok 

or simply {F} = [£(u;, £)]{u} where 

{F} = {Ni, F yv Vi, Mi] N 2 , F y2 , V 2 , M 2 } t , {u} = {u u v u w u u 2 , v 2 , w 2 , ip 2 } T 

and [ k ] is the [8 x 8] dynamic element stiffness matrix. This stiffness matrix is fre- 
quency and wavenumber dependent, complex, and symmetric. The stiffness relation 
for the throw-off or one-noded element is simply 

{F}i = = [£(o;,£)]{u}i (2.16) 

The stiffness matrix in this case is of size [4x4], These stiffness relations can be used 
to recover the results already presented in the cited references. 

For illustrative purposes, Figure 2.4 shows the normalized ku, k 22 , k 33 , and fc 44 
diagonal terms for the same shell used to illustrate the shape functions. As is typical 
with spectral elements, they exhibit a very large dynamic range. The normalizations 
are with respect to the stiffnesses for straight thin beams (Doyle, 1991) but modified 
for plates. That is, they are presented as 



kul{Eh/L ) , k 22 /(Gh/L), k 33 /(l2D/L 3 ) , fc 44 /(4D/L) 

Note that the kn stiffness is substantially less than the static values, and it is only 
after the cut-off frequency (~ 3 kHz for m = 0) does it become greater than unity. 
On the other hand, k 33 is always significantly larger than the static straight value; 
this is because the L 3 in the denominator predicts an inordinately small static value. 
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2.5 Point Excitation of a Shell 

The advantage of the curved element is that it can be combined with other elements 
— either flat or curved — to form significantly more complex structures. Danial and 
Doyle (1994, 1995) discuss aspects of the computer programming structure required 
for multiply connected spectral elements; they also show how the spectral element 
method can be hosted on a massively parallel machine. The variety of possibilities 
is too great to pursue here, so we will be content with two short examples. Both 
involve the point excitation of a complete cylindricjd shell but the latter shifts the 
boundaries so as to explore the nature of the reflections. 

We form the complete cylinder by combining two curved elements of the same 
material properties each of which is a half-circle. Because of symmetry, we could 
model the shell with a single half-circle element or as a single whole-circle element. 
In the latter case, half of the applied load must be placed at each node and the extra 
conditions of u\ = 0, = 0; U 2 — 0, ip 2 — 0 be imposed. The two element model 

was chosen so as to verify the assemblage process. 

The elements can be joined at the common nodes by merely summing together 
the appropriate dynamic stiffness matrix components to form a global stiffness matrix 
[ K ]. We must therefore first rotate each element stiffness to this global system. This 
must take into account that the ends of the curved element are oriented differently 
to each other and an appropriate scheme is illustrated by Bilodeau and Doyle (1997). 
The global stiffness matrix is then simply formed by adding the two global element 
stiffness matrices together. 

With the elements connected, we can now determine the responses due to a point 
excitation. The input force history used is the same as used by Danial et al (1996): it is 
a pulse of duration of about 120 /is, and has a frequency content of about 16 kHz. The 
velocity reconstructions for three points, 9 = 0°, 90", 180°, along the circumference 
are shown in Figure 2.5. For comparison purposes, the same shell was modeled using 
64 flat spectral elements. The two sets of results are indistinguishable even though 
there are very many reflections. It is worth pointing cut that when fewer flat elements 
were used the results deteriorated as the time increased. This shows the significant 
computational savings in using the curved element. 

There are obvious multiple reflections occurring The waves travel around the 
circumference and interact with each other. We can get an alternative insight by 
looking at the system response function G where u = GP and P is the input load. 
Note that this facility is an integral attribute of the spectral element formulation. 
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Prom Figure 2.6, we can see that many resonance frequencies are being established. 
These would be the frequencies of interest if an impulse/modal analysis experiment 
vibration experiment were being performed. It is interesting to note that the w 
response at 0° and 180° exhibit similar resonant behavior, but the response at 90° is 
missing many of the intermediate resonances. A significant peak appears at about 
3.3 kHz-, we now look further at this. 

We wish to explore the effect of the reflections on the formation of the spectral 
peaks. To do this, we will model the cylinder as a sequence of increasing elements but 
of the same radius — the cylinders can be viewed as forming helical coils. The net 
effect is to place Node 2 further and further away from Node 1. The final model is 
a single throw-off element. We use the same broad band excitation and the velocity 
reconstructions are shown in Figure 3.2. There are two points of interest. First, 
note how all traces have the same initial behavior — this is the duration before any 
reflections return. The 180° element shows the most reflections and these diminish 
with increasing element size. The second point of note the trailing oscillatory behavior 
even for the infinite element. Although the force excitation lasts less than 200 ps, the 
plate continues to oscillate in an almost resonant like fashion. Actually, a standing 
wave has been established. 

A different view of these behaviors is obtained by looking at the velocity responses 
in the frequency domain. It is clear in Figure 2.8 that there is the formation of 
an increasing number of spectral peaks as the element length is increased. This is 
expected but what is interesting is a peak in the ru-velocity response in the infinite 
length case. An analysis of the amplitude ratios shows that this peak is coming 
entirely from the first mode even though the response overall is dominated by the 
first mode. What is happening is that the curvature acts effectively as a continuous 
boundary and sets up a standing wave. 

The flexibility of the spectral element approach has allowed us to view the problem 
of the impact of a cylinder as one of wave propagation or vibration. This comes about 
because of the convenience in alternating between the time and frequency domains. 
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2.6 Discussion 

A spectral shell element was developed that extends the variety of problems the 
spectral element method can handle. It was shown to be accurate and certainly 
more computationally efficient than using multiple flat plate elements. The element 
developed possesses all of the features of the spectral element method which make the 
method desirable for solving dynamic problems. Primary among these are that the 
element can be very long and that the frequency domain formulation allows the system 
response functions to be determined automatically. This latter attribute, along with 
the fast Fourier transform, enables a duality between the time domain and frequency 
domain to be presented conveniently. Information not readily seen in one is often 
detectable in the other. The example problems demonstrate these features. 

One of the challenges that arose during the development was the problem of 
determining the spectrum relations. Even though formulas for the solution of quartic 
equations were derived, it is not always certain which branch should be chosen after a 
branch point occurs. In the quadratic case, which occurs for flat plates, the ambiguity 
is removed by adding damping. This becomes essential in the general case so that 
the phase of each kj(u>, £) can be tracked correctly — the effect of the damping is to 
separate each of the modes. The damping is also necessary in order to distinguish 
between the forward and backward moving waves. 

Based on the problems we have considered, it is clear that once there is more than 
one connection in a structure, it is essential to have a matrix methodology to handle 
the many unknowns. The spectral element approach presents itself as a well founded 
matrix method that embodies a number of efficiencies we have long associated with 
the conventional finite element method. For the range of problems they are suited 
for, they show great efficiencies and conveniences. Being formulated in the frequency 
domain means it is also ideally suited for energy flow analysis of the type described 
by Langley (1994). It is also appropriately formulated for tackling the solid/fluid 
interaction problems as occur in structural acoustics. A beginning in this direction is 
described by Bilodeau (1995). 

The approach as presented so far has difficulty vith localized discontinuities of 
properties or geometry unless the have a very simple geometry. At present, there 
is also the restriction that lateral boundaries must be relatively far away. Similarly, 
there is a restriction that the material properties and geometry in the lateral direction 
be uniform. It is hoped to tackle these issues in the future. 
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Figure 2.1: Coordinate system and displacements for a segment of cylindrical shell. 
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Figure 2.2: Spectrum relations kjh for a shell segment with R = 254 mm and h = 
25.4 mm, and m = 40. 
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Figure 2.5: Velocity reconstructions for the point excitation of a closed cylinder. 
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Figure 2.8: Frequency domain response at y = 2R for the point excitation of shells 
of different circumference. 





Chapter 3 

Acoustic Radiation from Plate 
Structures 


Spectral elements, which incorporate the effect of fluid loading on the structure, are 
developed for analyzing acoustic radiation from dynamically loaded extended plates. 
These elements may be conveniently joined to form complex thin-walled structures 
composed of many segments. 
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3.1 Introduction 

There axe many practical situations where the interaction between the dynamics of a 
structure and a surrounding fluid is of great importance. The most obvious is noise; 
noise is the propagation of acoustic energy through the fluid. The interaction can also 
influence the response of the structure itself; examples include dams, chimney stacks, 
ships, fuselages, propellers, and transmission cables. Fluid Loading problems are very 
hard to solve exactly, and for geometries and configurations of practical interest it 
is essential to be able to make useful simplifying approximations [31]. The purpose 
of this paper is to introduce a method being developed for analyzing folded plate 
structures immersed in a fluid. 

A schematic of the cross-section of the folded plate structures of interest is shown 
in Figure 3.1a; the plates extend in the y-direction. Such a structure when immersed 
in a fluid can experience three types of loading. The first is the structure-borne 
excitation caused by the propagating structural waves. The second is a pressure 
loading arising from some source within the fluid; this acts as a distributed external 
loading on the structure. The third load is the ‘self-loading’ or fluid loading caused 
by the moving structure interacting with the fluid. This also acts as a distributed 
pressure loading but, since the magnitude depends on the motion of the structure, it 
has the effect of coupling the structure and fluid motions. Additionally, these enclosed 
structures are viewed as having two distinct regions; an interior region and an exterior 
region. The interior region has loadings due to both the fluid and the radiation from 
the vibration of other plates, whereas, the exterior region only experiences loadings 
from the fluid. In this paper, we restrict the fluid to be only on the exterior of the 
structure so that we can concentrate on the structural dynamics, fluid loading, and 
radiation into the fluid. The reverberation problem f)r enclosed spaces is the subject 
of a later paper. 

Our analysis of the structural dynamics is based on the spectral element method 
(References [6, 8, 11 ] give summaries of the approach as applied to frame structures) 
but applied to structures of extended areas. This application to thin-walled struc- 
tures begins by combining the spectral analysis for i i-plane wave responses [19] and 
out-of-plane flexural behavior [27] to form a matrix method approach for folded plate 
structures [22], The variety of structural shapes enccmpassed was recently enhanced 
by adding a segmented curved shell element [32]. The method has some very use- 
ful characteristics. First, it is a frequency domain fcrmulation. Second, it is geared 
toward handling very large areas. That is, it is a counterpoint to conventional ele- 



56 


merits which always model a region as a collection of very small sub-domains. Third, 
it is stiffness formulated. This means the process of assemblage to form complex 
structures is relatively easy. 

The solid/fluid interaction problem is very intricate [33, 34, 35]; many of the finer 
points are covered in the excellent summary paper by Crighton [36]. The essential 
difficulty is that the two media are coupled in a convolution sense — this is unlike a 
plate on an elastic foundation, say. We tackle the interaction problem by incorporat- 
ing the effect of fluid loading on the structure directly into the element formulation. 
In this way, the structural formulation for in vacuum dynamic response is unaffected. 
Our approach is approximate but is shown to be reasonably accurate for medium 
fluid loading. Radiation from these extended plates is handled very conveniently 
by utilizing the shape functions associated with the spectral elements. In fact, this 
computation of the pressure response in the fluid is performed as a post-processing 
operation. The challenge we have here is to match the motion of the finite plate 
to that of the fluid. But the domain for the fluid is (at least) the half space above 
the plate and is considerably larger than the length of the plate itself. Therefore, to 
match the plate and fluid boundaries, we must extend the plate boundary. We do 
this by assuming the displacements can be matched by imposing w = 0 outside of 
the finite plate. We further assume that this can always be done even if the plate is 
not physically baffled but is attached to other plate segments. The idea is illustrated 
in Figure 3.1b. 

Our interest in this paper is on aspects of the structure/fluid interaction, therefore 
for simplicity, we take the plates and loading as being uniform in the ^/-direction. We 
also have in mind aeronautical structures as typified by aluminum in air — this corre- 
sponds to a relatively light fluid loading situation as compared to naval applications. 
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3.2 Plate/Fluid Interaction: Infinite Plate 

As a prelude to considering finite plates, we first look at the response of an infinite 
plate with fluid loading on one side. This will allow as to establish the context of the 
approximations to be made in the next section. The plate lies in the x — y plane, but 
for simplicity we take the plates and loading as being uniform in the y-direction. 

The governing equations for the deflection w(x, t) of a thin plate and the resulting 
pressure p(x,t), when a mechanical loading q(x, t) is applied, are 

d*w . dw , d 2 w . . . „ d 2 v 

D ~Q^4 + ' nh ~^ + ph ~d^' = q ^ x ’ t ^~ p ^ x ' z = 0 ’ t ^ P 


BV 2 p = p a - 


where h is the plate thickness, D the plate stiffness, p the plate mass density, B the 
fluid bulk modulus, and p a the fluid mass density. We have also added some damping 
to the structure through 77. The interface conditions between the plate and fluid are 


Wplate(x) = Wfl uid (x, 2 = 0 ), 


dp(x, 2=0) 

dz 


d 2 w(x, 2 = 0 ) 


Our solution technique uses spectral analysis; this assumes solutions of the form 


>(x, t) = Y w{x, u> n )e lWnt , p(x, z,t) = p(x, 2 , u) n )e u 


n2n . 

(3.3) 


where T is the time window for the discrete transform. Typically, N ranges as 
512 w 4096. When these are substituted into the governing differential equations, 
we get 

'll) { x 1 

D— — phu> 2 w{x) = q(x) - p(x, 2 = 0), BV 2 p + p a u) 2 p = 0 (3.4) 

where u 2 = u 2 — iut]/ p and the interface relations become 

t \ ■> 1 9p(x, 2 = 0) 9 . . . , 

Wplate\X) — Wfluid[X, Z — 0) , ~ — CO p a w{x ) (3-5) 

o z 

This set of equations will be our primary equations governing the dynamic response 
of the plate and fluid. 

Consider a single infinite sheet with fluid only on one side and with a line loading 
along x = 0. Let the plate deflection and loading be represented in the form 

M Mr, 

w(x) = Y, Wme~ lUx , q (z) = Y ^ » U = (3.6) 

m tti ** 

where W is the space window for the discrete trarsform. Typically, M ranges as 
500 « 2000. The pressure has a similar representation, except that we must realize 
it is two-dimensional 


P(x , z) = Y Pme iUx e ik ‘ z , k z = Jk* - 


(3.7) 
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This corresponds to pressure waves radiating from the plate surface in the positive 
z-direction. The pressure boundary condition at z = 0 is now 

-ikzPm - Pa^Wm 


Substituting these into Equation (3.4) gives the displacement of the plate and the 
pressure in the fluid as, respectively, 


w(x) 


p(x, z) 


= E 

m 

= E 


q m e * mX 

2 

i D tm ~ P huj2 + iuj r]h - ] 


ik 7 




-iZmX 


i D tm ~ P 1 ™ 2 + iuJ vh - ] ■ lkz 


ik z p a u) 2 


(3.8) 


As with the case of a plate in a vacuum, q m is chosen to be 1.0 to represent a line 
load at x = 0. 

Figure 3.2a shows the resulting pressures 100 mm from a 2.5 mm thick aluminum 
plate subjected to a line loading. The history of the loading is a smoothed triangular 
pulse of duration about 150 ps and having frequency content of about 16 kHz. The 
pressures exhibit an oscillatory behavior where the period of the zero crossings is 
almost constant. Figure 3.2b looks at the system response function |(j| where the 
pressure is related to the applied resultant load as p = GP. The most striking feature 
of the pressure is the spectral peak in the vicinity of 5 kHz — this corresponds to 
the coincidence frequency as discussed next. Note that the peak gets sharper further 
away from the impact site, and that there is more filtering before the coincidence 
frequency. 

The coincidence frequency occurs when the phase speeds of the plate bending 
wave and of the acoustic wave in the fluid are equal [33, 34] and is given by 


u) c 




yfi 2(1 - ^ 2 ) , 


Co 




(3.9) 


The significance of the coincidence frequency is that wave components in the plate 
above this frequency are radiated easily into the fluid as seen in Figure 3.2. The 
coincidence phenomena depends on the extent of fluid loading, thus it is useful to 
have a measure of this loading. Following Reference [36], fluid loading is charac- 
terized by two independent parameters: a mass ratio a = p a c a j phto and a speed 
ratio M = c/ c a = k a /k. We have M = 1 at coincidence. Both parameters are fre- 
quency dependent; we can arrange for only one parameter which is varied with uj by 
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introducing 

e = Ms. = EiEi/Ji2(i - * 2 ) 

phu c p c a v 

The parameter e, called the intrinsic fluid loading parameter, is the same for all plates 
of a given material embedded in a given fluid. This is typically small with values such 
as 

steel/water: 0.130, aluminum/air: 0.002 
In the examples that follow, we are primarily interested in the aluminum/air case. 
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3.3 Waveguide Modeling with Fluid Loading 


Our goal here is to formulate the structural dynamics problem in terms of a waveguide. 
If we can replace the wavenumber transform solution for the plate response given in 
Equation (3.8) with a waveguide solution then that saves a summation operation. 
But much more importantly, it allows us to terminate the waveguide and thereby we 
are in a position to assemble complex structures. The difficulty is that while the effect 
of the fluid loading is that of a distributed pressure, it’s non-local character makes 
it different from the pressure caused by a distributed spring, say, and hence we must 
invoke special procedures. 

To begin our construction of a plate waveguide, we ask if it is possible for free 
waves to propagate in the plate immersed in a fluid. That is, we seek wave solution for 
the plate of the form w = we~ lkx when the loading q is zero. This implies a pressure 
response of 


pe~ ikx e~ ikzZ 


\fkl - , 


k 2 

K a ~ ~2 

C a 


It must be borne in mind that as long as k z is chosen as above then irrespective of 
the value of k the fluid equations are satisfied. When both w and p are substituted 
into the governing equation we obtain the characteristic equation 


(k 2 - + 0 1 ) - ^ =0, H 2 =\J^ (3.10) 

The third term in the first equation describes the effect of the fluid loading. A detailed 
explanation of the significance of each term is given by Crighton [36] and the roots of 
this characteristic equation have been studied extensively in References [37, 38, 31]. 
Contour plots of the characteristic equation are shown in Figure 3.3. It is seen that 
there are four dominant roots similar to a plane in a plate. 

When the response of Equation (3.8) is viewed as a contour integral [39] in the 
complex planes of Figure 3.3, it has contributions from the poles and a contribution 
which arises from the branch cut associated with k z . This latter contribution is most 
significant at impact points and corners [36], By neglecting this contribution, we 
would then in a position to replace the responses with just the pole contributions. 
That is, we will have a waveguide representation. To quantify this contribution, we 
consider the response of a 25 mm thick steel plate in a fluid of density 138 g/m 3 and 
modulus 0.37 GPa subjected to the pulse line loading. This gives an intrinsic fluid 
loading of c = 0.020 which is an order of magnitude larger than the cases of actual 
interest. Figure 3.4 shows the responses with and without the branch cut contribution. 
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The agreement of the two solutions is very good with some deviations occurring at 
the load site that are primarily in the low frequency range. This is confirmed at the 
large x locations. We therefore conclude that in the case of an aluminum plate in air 
the branch cut contribution can be neglected. 

Numerical solutions for the roots of the characteristic equation are shown in Fig- 
ure 3.5 as circles along with the in vacuum roots indicated as the dashed lines. These 
plots, which are for 2.5 mm aluminum plates in air, show that the fluid loading pri- 
marily has the effect of altering the imaginary part of the first mode. We see that at 
coincidence and beyond, the effect is of increased viscous damping — this is consistent 
with the observation that the fluid is receiving more of the energy at these frequen- 
cies. But otherwise the behavior is very similar to the in vacuum behavior. That is, 
the pole contributions are associated with root ki which corresponds predominantly 
to the propagating flexural wave, and with root k 2 which corresponds predominantly 
to an evanescent flexural wave. 

Making the assumption that there are only two dominant structural waves allows 
us to obtain approximate analytical expressions for 1 he roots as 

1 = \[kl ~ I 2 

k, 2 = yfkl T(P (3. II) 

These approximations are also plotted in Figure 3.5. There is very good agreement 
with the exact numerical roots over the entire frequency range including the region 
near coincidence. 

By incorporating the fluid loading term directly into the modified spectrum rela- 
tions, we are now in a position to replace the double summation wavenumber trans- 
form solution by a single summation over frequency. This is useful when an enclosed 
structure is viewed as having two distinct regions; an interior region and an exterior 
region. The interior region has loadings due to both the fluid and the radiation from 
the vibration of other plates. However, the exterior egion only experiences loadings 
from the fluid. Using the modified spectrum relations, the fluid loading for the ex- 
terior problem can be accounted for without considering the fluid response. In this 
way, the solid/fluid interaction problem is partially decoupled. 

We conclude this section by illustrating the diffei ence the waveguide formulation 
makes. The transverse displacement of one half of th ; infinite plate with two forward 
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propagating waves is expressed as 

w(x) = Ae~ iklX + Be-* 21 


where ki and k 2 are the modified spectrum relations presented in Equation (3.11). 
The boundary conditions for this problem are that at x — 0 the slope of the plate is 
zero ( dw/dx = 0) and the applied load is related to the shear by 



V = -D 


cPw 

dx 3 


The response of the plate is then determined to be 


ie(x) 


P 

2Dik\ (k\ - kl) 


e ~ iklX 



(3.12) 


where most quantities are frequency dependence. This solution is to be compared to 
Equation (3.8). The responses of Figure 3.4 were computed using this solution. 

It is clear that the waveguide approximate modeling for plates in a surrounding 
fluid is accurate and efficient. 
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3.4 Modeling of Finite Plates in a Fluid 

We now illustrate how the waveguide modeling allows us to tackle the very difficult 
problem of the response of finite plates in a fluid. 

We begin by writing the general transverse displacement for a typical plate seg- 
ment of length L in the form 

w(x) = Ae~ iklX + Be~ ik2X + Ce~ ikl(L ~ x) + De -<fca(L-x) 


where A, B, C and D are constants to be determined from the boundary conditions 
on the segment. Also, by using the modified spectrum relations, we claim that this 
adequately represents the behavior of a plate immersed in a fluid. It is advantageous, 
when dealing with finite or multiply connected structures, to use a solution formu- 
lation that already incorporates the connectivities. The end conditions on the plate 
segment are 


w(0) = W\ , 


dw( 0) 


= 0i , u>(L) = w 2 , 


dw(L) 


02 


dx ’ dx 

Solving for the coefficients in terms of the nodal deg ees of freedom allows the trans- 
verse displacement of the plate to be re-written in the form of a collection of shape 
functions 


w(x) = 9 i{x)wi+ g 2 (x)LTpi+ g 3 (x)w 2 + g 4 (x)Lfo (3.13) 

The frequency dependent shape functions gj(x,u) are given as 

9\{x) = [r\'h\{x) + r 2 h 2 (x )]/ A 

g 2 {x) = [rih 3 {x) + r 2 hi(x)]/ A 

<7a(z) = [rMx) + r 2 h 1 {x)]/A (3.14) 

g A {x) = -\r x h A (x) +r 2 h 3 (x)]/ A 

t i = i(k x — Ar 2 )[l — e~ klL e~ k2L ] , r 2 = i(k i + A: 2 )[l — e~ klL e~ k2L ] 

where A = —(r\ + r%) and 

h\(x ) = +ik 2 [e~ lklX — e~ lk2L e~ ,kl ^ L ~ x )} — iki\e~ lk2X — e~ iklL e~ ik2 ^ L ~ x ^] 

h 2 (x) = —ik 2 [e~ lk2L e~ lklX — -j. ik x \e~ lklL e~ ik2X — 

h 3 (x) = [ e -'b* + e -it 2 £ e -i*i(t-i)j _ \ e ~ik 2 x + e -ikiL e -ik 2 (L-x)-\ 

h\{x) — ^ e ~ik 2 L e -ikix e ~iki{L-x)^ _ jg-tkiLg-ikjz e -ik 2 {L-x)^ 

As an example, the shape functions g x and g 2 are shown plotted in Figure 3.6 for 
a number of frequencies. The shape functions occur in pairs where g 3 and g A are 
the mirror images of g x and g 2 , respectively. These shape functions are comparable 
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to those of the conventional finite element method except that they are frequency 
dependent and can represent very large areas. References [11, 32] illustrate how 
these are used as the basis for the spectral element representation of the folded plate 
structures, these references also show the assemblage procedures. 

What we want to look at here is how a typical plate segment radiates pressure 
into the fluid. The challenge we have here is to match the motion of the finite plate 
to that of the fluid. But the domain for the fluid is (at least) the half space z > 0 and 
— oo < x < oo which is considerably larger than the length of the plate. Therefore, 
to match the plate and fluid boundaries, we must extend the plate boundary in the 
x direction. If the finite plate is baffled, that is, extended on both sides with very 
stiff material, then the displacements can be matched by imposing w = 0 outside of 
the finite plate. We will assume that this can always be done even if the plate is not 
physically baffled but is attached to other plate segments. The schematic is shown in 
Figure 3.1b. 

At the surface of the plate, 2 = 0, the fluid displacement must be equal to the 
plate displacement. By extending the plate deflection over the full space window of 
the fluid, we can then give it the spectral representation 


w(x) = 

m 

Applying this to the shape functions gives 


= Widim + Ltj)ig 2m + w 2 gz m + Ltp2g im , 9jm= I gj{x)e lUx dx (3.15) 

J w 

These integrals are easily evaluated since gj{x) contain only exponentials. The shape 
functions gj(x) are zero outside the length of the plate element; this is equivalent to 
assuming each finite plate segment is baffled to infinity. The continuity of Wj and ipj 
between segments ensure continuity of the fluid field. The response of the fluid now 
has the spectral representation 


w{x, z) = 53 


w r 


e ~ ik z z e -itmX 


p(x,z) 


Pa W 

m 


—ik : 


w m e~ ikzZ e~ 


i£mX 




These, in combination with Equation (3.15), relate the fluid response at any point to 
the plate nodal degrees of freedom. If there are many plate segments, then the total 
response would be the sum of the contributions from each segment. We will illustrate 
these formulas with the example of a finite plate. 

In piecing together the solution for the finite plate of Figure 3.1c, we treat the 
plate as having two segments with the boundary conditions 


segment 12: 


w x = 0 ip i = 0 ; w 2 = w a , = 0 
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segment 23: w\ = w a , tj)\ = 0; w 2 = 0ip 2 = 0 

This leads to the solution for the fluid displacement written in terms of the central 
deflection of the plate as 

w(x, z) = u> 0 £ [hme-*^ x+L) + gi m e-* mX ] e~ ikzZ 

m 

Note that the representation for the first plate segment must be shifted an amount 
L. The pressure is given by 

p(x, z )=WoY J ^ Le-^ (l+i > + e~ ikzZ 

m l ^z 

The sheax relation at the load location is \P — —V = Dw 0 "'. This gives the central 
deflection as 

, _ P 

W ° ~ 2£>p 1 "'(0) 

With reference to Figure 3.1c, it is worth emphasizing that although w is physically 
continuous at Node 2, the above representations have a discontinuity. That is, seg- 
ment 12 is discontinuous to the right while Segment 23 is discontinuous to the left. 
Hence this example is a good first test case of our scheme. 

To test the validity of our approximations, a planar finite element model of a 
baffled plate in a fluid was constructed for the problem shown in Figure 3.1c; only 
a finite fluid domain was modeled. The plate is 25 mm thick steel and the fluid has 
the properties such that e = 0.02. Figure 3.7 shows a comparison of velocities of the 
plate and the fluid. The agreement is quite good up to the time when reflections come 
from the far boundaries of the fluid — the spectral solution has no such boundaries. 
The agreement is especially good for our purpose when it is realized that the fluid 
loading corresponds to a loading factor nearly ten times that of our cases of interest 
of aluminum plates in air. 

These equations were also applied to the case of a finite 2.5 mm thick aluminum 
plate in air. The responses for the plate are showr in Figure 3.8; the presence of 
multiple reflections are obvious. Figure 3.9 shows the corresponding frequency domain 
behavior where the reflections give rise to multiple spectral peaks. Also shown are 
the resonance frequencies for a vibrating plate in vacuum, it is clear that the impact 
has excited many of the symmetric modes of vibration. The pressure responses for 
the fluid can also be seen in Figure 3.8. Similar to ti e plate response, these indicate 
the presence of multiple reflections occurring in the date. This figure also indicates 
that the baffled finite plate not only excites the fluid directly in front of it but also at 
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a distance along the baffle. The pressure at x — 2 L, z = 0 is non-zero even though 
the plate is baffled at that location — this is further indication that, in fluids, the 
relation between the pressure and displacement is non-local. The frequency domain 
depiction of these responses is shown in Figure 3.9. It can be seen that the frequencies 
at which the structural resonances were present in Figure 3.9 are readily transmitted 
into the fluid. 

Intuitively, we might have thought that the response should be largest along a 
line normal to the plate. The responses shown in Figure 3.8 indicate that this is not 
so. Furthermore, Figure 3.9 shows that it is even frequency dependent. Figure 3.10 
shows the directivity patterns at a number of frequencies; the near field behavior was 
computed from the full solution with r = 2 L. It is clear that these patterns are very 
sensitive to direction when the frequency gets close to coincidence. Also shown are 
the deflected shapes of the plate at each frequency. These shapes indicate an almost 
sinusoidal plate deflection except at the center and edges — these are the points of 
significant radiation. 
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3.5 Discussion 

The methods discussed in this paper is a first step in extending the spectral element 
method to include structure/fluid interaction problems. The key step is that by 
incorporating the fluid loading into the spectrum relations allows us to maintain 
the element formulation. This is important because the spectral element approach 
presents itself as a well founded matrix method that embodies a number of efficiencies 
we have long associated with the conventional finite element method. For the range of 
problems they are suited for, the spectral elements have been shown to conveniently 
model wave propagation in structures made of multiple panels. Furthermore, the 
spectral element is eminently suited for hosting on massively parallel computers [28, 
21 ]- 

There are many more developments needed. The most important of them are: 
Verify the radiation approximation from re-entrant edges, Implement radiation from 
curved surfaces, and, Implement the book-keeping necessary for solving reverberation 
in a enclosed space. We will to report on these developments in future papers. 
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Figure 3.1: Geometries of the plate structures with fluid loading, (a) Schematic of 
cross-section of folded plate structures of interest, (b) Treatment of plate connections 
as baffled extensions, (c) Modeling of a finite baffled plate as two plate segments. 
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Figure 3.2: Time and frequency domain responses for the line loading of an infinite 
plate, (a) Fluid pressures p(x, z = L,t). (b) System response function \G(x,z = L)\ 
for pressure. 
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Figure 3.3: Contours of the characteristic equation for two frequencies straddling 
coincidence. 
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Figure 3.4: Comparison of the wavenumber trans orm solution to the waveguide 
solution that neglects the branch cut contribution. 



72 



Figure 3.5: Wavenumber behavior for an aluminum plate in air. Left is ki(u), right 
is k 2 (u>). 
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Figure 3.6: Shape functions < 7 i(x) (top) and g'^x) (bottom) for a number of frequen- 
cies. 
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Figure 3.7: Response with modified spectrum relations and comparison with finite 
element solution. 
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Figure 3.8: Time domain responses for a line loaded aluminum plate of length 2 L = 
500 mm. (a) Plate responses w(x,t), (b) fluid pressures p(x, z = L,t). 
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Figure 3.9: Frequency domain responses for an impacted plate of length 2 L = 
500mm. (a) Plate responses \iuw(x)\, (b) fluid pressure | p(x,z — L ) |. 
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Figure 3.10: Directivity patterns for impacted finite plate. Also shown are the plate 
shapes at each frequency. 
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